Rcpp::sourceCpp("/home/fra/Documents/GitHub/CommonAtomModel/Functions/CAM/CAM.cpp")
source("/home/fra/Documents/GitHub/CommonAtomModel/Functions/CAM/newCAM.R")
## 
## Attaching package: 'reshape'
## The following objects are masked from 'package:plyr':
## 
##     rename, round_any
## Registered S3 method overwritten by 'pryr':
##   method      from
##   print.bytes Rcpp
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.4     ✓ dplyr   1.0.2
## ✓ tidyr   1.1.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::arrange()   masks plyr::arrange()
## x dplyr::combine()   masks gridExtra::combine()
## x purrr::compact()   masks plyr::compact()
## x purrr::compose()   masks pryr::compose()
## x dplyr::count()     masks plyr::count()
## x tidyr::expand()    masks reshape::expand()
## x dplyr::failwith()  masks plyr::failwith()
## x dplyr::filter()    masks stats::filter()
## x dplyr::id()        masks plyr::id()
## x dplyr::lag()       masks stats::lag()
## x dplyr::mutate()    masks plyr::mutate()
## x purrr::partial()   masks pryr::partial()
## x dplyr::rename()    masks reshape::rename(), plyr::rename()
## x dplyr::summarise() masks plyr::summarise()
## x dplyr::summarize() masks plyr::summarize()
## Loading required package: coda
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## ##
## ## Markov Chain Monte Carlo Package (MCMCpack)
## ## Copyright (C) 2003-2021 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park
## ##
## ## Support provided by the U.S. National Science Foundation
## ## (Grants SES-0350646 and SES-0350613)
## ##
## 
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
## 
##     discard
## The following object is masked from 'package:readr':
## 
##     col_factor

Data

6 distribution osservate, 4 gruppi

x0 <- replicate(1, c(rnorm(50,3),rnorm( 150,10,.5)))
x1 <- replicate(1, c(rnorm(100,3),rnorm( 100,10,.5)))
x2 <- replicate(3, c(rnorm(100,-5),rnorm(100,10)))
x3 <- c(rnorm(100,3),rnorm( 100,10,.5),rnorm( 100,-5,.5),rnorm( 100,0,.5))

yd <- c(x0,x1,x2,x3)
yg <- rep(1:6,c(rep(200,5),400))
  
plot(yd,col=yg)

Run the model

m1 <- CAM(y_obser = yd,
            y_group = yg,
            K0 = 20,
            L0 = 20,
            prior = list(
              # hyperparameters NIG
              m0= mean(yd), 
              k0=1/var(yd), 
              a0=3, b0=1,
              # hyperparameters alpha and beta
              a_alpha=3, b_alpha =3,
              a_beta =3, b_beta = 3),
            nsim = 20000,
            burn_in = 20000,
            thinning = 1,
            verbose = F,
            fixedAB = F,
            kappa=0.5,
            cheap=F,
            seed=123456,sigma.prop.beta = 1
  )

Alpha e Beta

plot(ts(m1$A_DP)) 

acf(m1$A_DP)

pacf(m1$A_DP)

plot(ts(m1$B_DP))  

acf(m1$B_DP)

pacf(m1$B_DP)

Clustering

r1=mcclust::comp.psm(as.matrix(t(map_dfc(m1$Z_j, ~.x)  )))
## New names:
## * NA -> ...1
## * NA -> ...2
## * NA -> ...3
## * NA -> ...4
## * NA -> ...5
## * ...
image(r1)

a1=mcclust.ext::minVI(r1,method = "greedy")
plot(yd,col=rep(a1$cl,table(yg)))

r=mcclust::comp.psm(as.matrix(t(map_dfc(m1$Csi_ij, ~.x)  )))
## New names:
## * NA -> ...1
## * NA -> ...2
## * NA -> ...3
## * NA -> ...4
## * NA -> ...5
## * ...
image(r)

a=mcclust.ext::minVI(r)


plot(yd,col=a$cl)

Densities

xx <- seq(-10,15,by = .05)
LLL <- list()

LLL <- posterior.densities(m1,howfarback = 1000,T,xx = xx)
## 123456
par(mfrow=c(1,2))
for(hh in 1:6){
hist(yd[yg==hh],breaks = 30,freq=F, main="hist+mcmc+mean")
for(i in 1:1000)
lines(LLL[[hh]][i,]~xx,col=1)
lines(colMeans(LLL[[hh]])~xx,col=2,lwd=3)

hist(yd[yg==hh],breaks = 30,freq=F, main="hist+mean")
lines(colMeans(LLL[[hh]])~xx,col=2,lwd=3)
}

# sim <- 500
# ind <- 500
# ygind <- yg[[ind]]
# 
# plot(yd)
# points(yd[ind]~c(ind),col=2,pch=21,bg=2)
# 
# mat <- matrix(NA,nsim,2)
# for(sim in 2:nsim){
# 
# mat[sim,1] <- m1$OMEGA[[sim]][ m1$Csi_ij[[sim]][ind] 
#                  ,
#                   m1$Z_j[[sim]]  [ygind]]
# mat[sim,2] <- m1$THETA[[sim]][ m1$Csi_ij[[sim]][ind] 
#                  ,
#                  1]
# }
# 
# plot(ts(mat))